R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter. # Tue Jun 7 19:35:18 2022 ——————————

source("tianfengRwrappers.R")

rat Frontier GSE174098 carotid

rat10x <- CreateSeuratObject(Read10X("./rat_scRNAseq/"), names.field = 2, names.delim = "-",
                                   project = "rat", min.cells = 10, min.features = 300) %>%
  PercentageFeatureSet(pattern = "^Mt-", col.name = "percent.mt") 

table(rat10x$orig.ident)
VlnPlot(rat10x,"nCount_RNA") /
VlnPlot(rat10x,"percent.mt") /
VlnPlot(rat10x, "nFeature_RNA")


rat10x <- rat10x %>% subset(subset = nFeature_RNA > 400 & nFeature_RNA < 4000 &
                              nCount_RNA > 1000 &  nCount_RNA < 30000 & percent.mt< 10) %>%
    SCTransform(vars.to.regress = "percent.mt", verbose = F) %>% 
    RunPCA() %>% FindNeighbors(dims = 1:20) %>% 
    RunUMAP(dims = 1:20) %>% 
    FindClusters(resolution = 0.1)

mouse coronary GSE131778

mouse_coronary_countmatrix <- read.csv("./GSE131776_mouse_scRNAseq.txt", sep = "\t")
func <- function(s) {
  paste0(strsplit(s, ".", fixed = T)[[1]][2], "_", strsplit(s, ".", fixed = T)[[1]][1])
}
colnames(mouse_coronary_countmatrix) <- lapply(colnames(mouse_coronary_countmatrix), func) # 拆分样本

mousecor <- CreateSeuratObject(counts = mouse_coronary_countmatrix,
                                   project = "mouse_cor", min.cells = 10, min.features = 300) %>%  PercentageFeatureSet(pattern = "^mt-", col.name = "percent.mt") 

# saveRDS(mousecor,"mousecor.rds")
table(mousecor$orig.ident)
VlnPlot(mousecor,"nCount_RNA") /
VlnPlot(mousecor,"percent.mt") /
VlnPlot(mousecor, "nFeature_RNA")


mousecor <- mousecor %>% subset(subset = nFeature_RNA > 400 & nFeature_RNA < 4000 &
                              nCount_RNA > 1000 &  nCount_RNA < 30000 & percent.mt < 10) %>%
    SCTransform(vars.to.regress = "percent.mt", verbose = F) %>% 
    RunPCA() %>% FindNeighbors(dims = 1:20) %>% 
    RunUMAP(dims = 1:20) %>% 
    FindClusters(resolution = 0.1)
table(mousecor$orig.ident)
saveRDS(mousecor,"mousecor.rds")

SMC2

mousecor <- readRDS("mousecor.rds")
f("Prdm16",mousecor, label.size = 7) + theme(legend.text = element_text(size = 20))

umapplot(mousecor_stromal, label.size = 6) + theme(legend.text = element_text(size = 20))

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the
existing scale.

mouse carotid scRNAseq GSE155513

dataload

count_mats <- list.files("./GSE155513_RAW/")
count_mats <- count_mats[count_mats != "sampleinfo.txt"]
allList <- lapply(count_mats, function(file) {
  dd <- read.table(paste0("./GSE155513_RAW/", file), row.names = 1,stringsAsFactors = F)
  colnames(dd) <- as.character(dd['gene',])
  dd <- dd[-1,]
  CreateSeuratObject(
    counts = dd,
    project = file, min.cells = 10, min.features = 300
  )
})
# 合并seurat对象
mouse_carotid <- merge(allList[[1]], 
  y = allList[-1], add.cell.ids = count_mats,
  project = "mouse_carotid"
)
rm(allList)

# saveRDS(mouse_carotid,"mouse_carotid.rds")
# View(mouse_carotid)

process

mouse_carotid <- readRDS("mouse_carotid.rds")

mouse_carotid <- mouse_carotid %>%  PercentageFeatureSet(pattern = "^mt-", col.name = "percent.mt") %>%
subset(subset = nFeature_RNA > 500 & nFeature_RNA < 3000 &
                              nCount_RNA > 1000 &  nCount_RNA < 20000 & percent.mt < 10) %>%
    SCTransform(vars.to.regress = "percent.mt", verbose = F) %>% 
    RunPCA() %>% FindNeighbors(dims = 1:20) %>% 
    RunUMAP(dims = 1:20) %>% 
    FindClusters(resolution = 0.1)

saveRDS(mouse_carotid,"mouse_carotid.rds")

SMC

stromal cells

mouse_carotid_stromal <- subset(mouse_carotid,idents = c(0,2,1,7))
mouse_carotid_stromal <- mouse_carotid_stromal %>% RunPCA() %>% FindNeighbors(dims = 1:20) %>% 
    RunUMAP(dims = 1:20) %>% FindClusters(resolution = 0.1)
mouse_carotid_stromal <- mouse_carotid_stromal %>% FindClusters(resolution = 0.2)

mouse_carotid_stromal <- readRDS("mouse_carotid_stromal.rds")
# saveRDS(mouse_carotid_stromal,"mouse_carotid_stromal.rds")

umapplot(mouse_carotid_stromal)
f("Dlx2",mouse_carotid_stromal) #Dlx2,Dlx5,Dlx6共同定位

分群表

table(group_tab[Dlx5poscells])/table(group_tab) 

GSM4705592_RPS003_matrix.txt.gz GSM4705593_RPS004_matrix.txt.gz GSM4705594_RPS011_matrix.txt.gz GSM4705595_RPS012_matrix.txt.gz GSM4705596_RPS007_matrix.txt.gz 
                    0.080645161                     0.076949502                     0.130000000                     0.029882604                     0.058282209 
GSM4705597_RPS008_matrix.txt.gz GSM4705598_RPS001_matrix.txt.gz GSM4705599_RPS002_matrix.txt.gz GSM4705600_RPS017_matrix.txt.gz GSM4705601_RPS018_matrix.txt.gz 
                    0.004746835                     0.039230575                     0.003976143                     0.093492209                     0.095634096 
GSM4705602_RPS013_matrix.txt.gz GSM4705603_RPS014_matrix.txt.gz GSM4705604_RPS015_matrix.txt.gz GSM4705605_RPS016_matrix.txt.gz 
                    0.066773504                     0.019255456                     0.042503503                     0.030303030 

图

human ds2

p <- multi_featureplot(c("FRZB","SOST","DLX5","DLX6"), ds2, labels = NA, label.size = 6)
ggsave("refds2_SMC2_carotid_stromal.png",plot = p, height = 7, width = 7,device = png)

Warning: Removed 1 rows containing missing values (geom_text).

human bulk RNA-seq GSE120521 carotid stable/unstable FPKM

fpkm2tpm <- function(fpkm){
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}

fpkm_matrix <- read.csv("GSE120521_FPKM.csv")
# fpkm_matrix <- distinct(fpkm_matrix) #去除重复行
fpkm_matrix <- fpkm_matrix[!duplicated(fpkm_matrix$name),]
rownames(fpkm_matrix) <- fpkm_matrix$name
fpkm_matrix$name <- NULL

tpm_matrix <- apply(fpkm_matrix, 2, fpkm2tpm)
colSums(tpm_matrix)

group_file <- c("stable","unstable","stable","unstable",
                "stable","unstable","stable","unstable")
boxplot(tpm_matrix, las = 2)

expr_mat <- tpm_matrix[!apply(tpm_matrix, 1, function(x){sum(floor(x) == 0)>3}),]

boxplot(expr_mat, las = 2)

library(limma)
expr_mat <- normalizeBetweenArrays(expr_mat)
expr_mat <- log2(expr_mat+1) #使用log2 scale

#PCA
library(ggfortify) 
df <- as.data.frame(t(expr_mat)) 
df$group <- group_file 
autoplot(prcomp(df[,1:(ncol(df)-1)]), data=df, colour = 'group')+ theme_bw() 

# If the sequencing depth is reasonably consistent across the RNA samples, then the simplest and most robust approach to differential exis to use limma-trend.

fit <- lmFit(expr_mat, group_file)
fit <- treat(fit, lfc=log2(1.2), trend=TRUE)
topTreat(fit, coef=ncol(design))

library(ggpubr)
dat <- expr_mat
design <- model.matrix(~factor(group_file))
fit <- lmFit(dat, design)
fit <- eBayes(fit)
# options(digits = 4)
topTable(fit,coef=2,adjust='BH')
deg <- topTable(fit,coef=2,adjust='BH',number = Inf)
head(deg) 

write.csv(deg,"./datatable/stable vs unstable.csv")
## look up FRZB, SOST, PRDM6, OGN
ds1markers[ds1markers$cluster == "SMC2",]$gene

ds2markers[ds2markers$cluster == "SMC3",]$gene

smc2markers <- intersect(ds1markers[ds1markers$cluster == "SMC2",]$gene, ds2markers[ds2markers$cluster == "SMC3",]$gene)


deg[intersect(smc2markers, rownames(deg)),]

ds1markers[ds1markers$cluster == "SMC1",]$gene

ds2markers[ds2markers$cluster == "SMC1",]$gene

SMC1markers <- intersect(ds1markers[ds1markers$cluster == "SMC1",]$gene, ds2markers[ds2markers$cluster == "SMC1",]$gene)

deg[intersect(SMC1markers, rownames(deg)),]

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.